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Abstract 


This work presents a solution to compute pressure change during injectivity tests in single-layer reservoirs 
with vertical wells using impulse functions (or Green’s functions). The proposed formulation employs a 
radially composite reservoir approach to depict the waterfront propagation and assumes a piston-like fluid 
displacement. Then, pressure change in multilayer commingled reservoirs may be obtained by combining 
the single-layer solution for each layer. Layer properties such as permeability, thickness and porosity may 
be different in each layer. The presented single-layer solution is also combined with Newmann’s product 
to derive an approximate formulation for pressure change in single-layer reservoirs with horizontal wells. 
Pressure and layer flow-rate profiles computed using the proposed formulation showed a good agreement 
when compared to previously existing solutions and a commercial flow simulator. Moreover, this work 
also estimates layer individual permeabilities from pressure and layer flow-rate data. 


Keywords: Injectivity Tests, Pressure Transient Analysis, Green’s Functions, Reservoir Characterization 





dealing with fluid disposal or flaring - hence, pre- 
senting lower environmental impact than conven- 


tional well testing (Borello et al.} |2019| 


1. Introduction 


Pressure transient analysis is crucial in reservoir 
engineering, as it provides valuable estimates about 


reservoir features and parameters. 
Typically, well testing is performed to estimate 


reservoir properties (Lefkovits et al.| 
1989} |Ahmed} |2010). Alternatively, parame- 


ters such as permeability and outer boundary con- 
dition may be estimated based on data collected 


during an injectivity test (Banerjee et al.| /1998| 
Peres et al, [2006). 


Injectivity tests have gained more attention re- 
cently, as they are operationally safer and avoid 
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2020). The reduced carbon footprint com- 


pared to conventional well testing becomes partic- 
ularly important considering that many oil field 
companies committed to substantially decrease 
greenhouse gas emissions in a foreseeable future 


Pressure transient data may be easily inter- 
preted using analytical solutions that accurately 
describe pressure behavior in the reservoir. Addi- 
tionally, analytical models are effortlessly coupled 
with nonlinear regression techniques and optimiza- 


tion methods (Coutinho et al} {2010} [Silva etal] 


Thompson and Reynolds} (1997) developed a 
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solution for single-layered radially heterogeneous 
reservoirs. ‘Their model has inspired the devel- 
opment of analytical solutions for pressure change 
during injectivity tests 
and Reynolds} [2009) [Bonafé et al [2020). 

(1998) presented a formulation 
for injection/falloff tests in single-layer reservoirs 
with vertical wells. As outlined by 
et al. (1998), their model assume that the rate- 
transient pulse propagates faster than the flood 
front. Hence, flow-rate inside the flooded region 
remains constant during the injection period, as 
illustrated in Fig. 

established objective 
criteria to verify that, in single-layer reservoirs, 
such condition holds for most cases of practical in- 
terest and also proposed a solution for injectivity 
tests in single-layer reservoirs with horizontal wells. 

inchuded. the 
thermal effects into the formulation and applied a 
non-linear regression method to estimate reservoir 
properties. provided an ap- 
proximate model for injectivity tests in single-layer 
reservoirs considering a multiple flow-rate scheme. 

Analytical solutions for injectivity tests in mul- 


tilayer reservoirs (Barreto Jr. et al.| 
2020b) were inspired by the single- 
layer solutions (Banerjee et al., |1998} 
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Figure 1: Schematics of rate transient front and water sat- 
uration front in a single-layer reservoir 












,|Peres et al.| , |Boughara and 
,|Bonafé et al.|!2020) and, hence, by 
the [Thompson and Reynolds} (1997) steady-state 


theory. However, in multilayer reservoirs, layer 
flow-rates at the wellbore may change in time due 


to differences between layer properties 1987, 
Ehlig-Economides and Josephj |1987| |Bela et al.| 
2020a |Silva et al.| 2021). 


Fig. |2} exhibits an illustrative representation of 
layer flow-rate profiles at the wellbore and the re- 
spective rate transient fronts in a two-layered reser- 
voir. As may be observed in Fig. |2| it is not possi- 
ble to assure that flow-rate is constant throughout 
the flooded region in multilayer systems, due to the 
rate transient that occurs at the wellbore. There- 
fore, in a strict sense, the [Thompson and Reynolds} 
steady-state theory does not hold for mul- 
tilayer formations. 

proposed an alternative solu- 
tion for injectivity tests in single-layer reservoirs. 
The flooded zone and the uninvaded region are 
described using a radially composite reservoir ap- 
proach, where the interface between regions (that 
is, the water front radius) moves in time 
and Horne} |1990). The formulation presented by 
Neto et al.| (2020) does not require flow-rate within 
the flooded region to be constant. On the other 
hand, it assumes a piston-like fluid displacement. 

also employed a radially 
heterogeneous reservoir model to develop a solu- 
tion for multilayer stratified reservoirs, while|Viana| 
included formation crossflow effects into the 
solution. 

Green’s functions (or impulse functions) for 
single-phase flow of slightly compressible fluids in 
single-layer reservoirs are known (Gringarten and) 
Rane 
(1991) developed an algorithm to deter- 
mine the pressure response in commingled reser- 
voirs using the appropriate single-layer impulse 
function for each layer. This algorithm was en- 


dorsed by (1994), in a work that also 


detailed how layer flow-rates at the wellbore may 
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Figure 2: Layer flow-rate profiles at the wellbore (on the left) and rate and saturation fronts in each layer 


be computed. (2019) applied Green’s 


functions to develop a formulation for well tests 
in two-layer reservoirs with formation crossflow. 
However, to the best of our knowledge, impulse 
functions have not yet been applied to evaluate the 
pressure change in reservoirs under two-phase flow. 


Therefore, the first goal of this work is to ap- 
ply Green’s functions to determine the pressure 
change during injectivity tests in single-layer reser- 
voirs using a radially composite approach. Then, 
the single-layer solution is coupled to the algorithm 


presented by |Kuchuk and Wilkinson) (1991) and 
Spath et al.| (1994) to obtain pressure and flow- 


rate profiles in multilayer reservoirs. Additionally, 
an approximate solution for pressure change during 
injectivity tests in single-layer reservoirs with hori- 


zontal wells is derived by applying /Newman| (1936) 


product and the presented single-layer formulation. 


The proposed models were validated via com- 
parison, for a set of synthetic cases, to a com- 
mercial finite difference-based flow simulator and 
to the previously existing solutions (Peres and| 


Reynolds} |2003}|Barreto Jr. et al.|/2011| 
Furthermore, the Delta Transient 


Method ((Galvao and Guimaraes} |2017| 
and Galvao\ ) was applied to estimate layer 
permeabilities in the multilayer cases. 








2. Proposed Formulation 


This Section outlines the proposed formulation 
for injectivity tests in single-layer reservoirs with 
vertical wells using impulse functions. This formu- 
lation may be coupled to the algorithm introduced 
by to evaluate the 
pressure behavior in multilayer reservoirs with ver- 
tical wells. Besides, an approximate model for in- 
jectivity tests in single-layer reservoirs with hor- 
izontal wells is obtained, combining the proposed 
model for radially composite reservoirs with an un- 
orthodox application of product. 

The reservoir is assumed to be homogeneous, lat- 
erally infinite and initially in equilibrium. Fluids 
are considered to be slightly compressible, with 
constant viscosity and compressibility. The pre- 
sented formulation considers that capillary and 
gravitational forces are negligible. An isothermal 
flow is considered and a piston-like fluid displace- 
ment is assumed. Finally, all computations are 
made assuming a consistent set of units. 


2.1. Impulse Function for Single-Layer Reservoirs 
with Vertical Wells 


Using a radially composite approach, water in- 
jection in a single-layer reservoir may be under- 
stood as a problem of a radially heterogeneous 
reservoir with a moving boundary rp(t) (Bratvold 


tion damage is included into the model by defin- 
ing an annular region, concentric to the wellbore 
axis, where the permeability k, is different than the 
reservoir permeability k 
[2020). Wellbore storage will be initially neglected 
and incorporated into the model afterwards. 

Figs. [3] to [6] show the reservoir top and lateral 
view at two distinct stages: while the water front is 
within the damaged region (Figs. [3] and [4) and af- 
ter the flooded region has overcome the skin radius 


(Viana) 2021) [Mastbanm et a, 202i): 


1 
/ 

Mil (2) 
where S,,; represents the water initial saturation 
and S,, denotes the oil residual saturation. 

The reservoir, then, may be split into three dis- 
tinct regions, depending on the flood front position. 
While the water front is within the damaged zone: 


Reg. 1: goes from r =r, to r = rp(t) 


(Figs. |5) and (6). In Figs. |3] to (6 the skin radius Reg. 2: goes from r = rp(t) tor =r, 


is represented as r,, while r,, denotes the wellbore 


radius, re denotes the water front radius and r, Reg. 3: goes from r =r, to re 


denotes the reservoir outer radius. 
The water front radius at each time step is de- 


termined as (Buckley and Leverett} |1942): 


j q(T )dr 


ie ah Fu lSw), (1) 





where qg represents the injection flow-rate, @ de- 
notes the reservoir porosity, h indicates the for- 
mation thickness and f’,(S,,) stands for the frac- 
tional flow derivative. For a piston-like fluid dis- 


placement, f/, is computed as (Neto et al. |2020| 





Figure 3: Reservoir top view for rp(t) < rs. 


And after the water front has overcome the skin 
zone: 


Reg. 1: goes fromr=ry tor=T; 
Reg. 2: goes from r =r, tor =rp(t) 


Reg. 3: goes from r = rp(t) to re 


According to the model assumptions previously 
mentioned, flow in each region obeys the diffusiv- 


ity equation (Ehlig-Economides and Joseph} |1987| 
2019): 
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Figure 4: Reservoir lateral view for rp(t) < rz. 





Figure 5: Reservoir top view for rp(t) > rs. 





r Or . Or ™, Ot 
where the subscript 7 € {1,2,3} indexes the re- 
gion, r stands for the radius, t denotes the time and 
Ap(r,t) = p(r,t)—p(r, t = 0) indicates the pressure 
change resulting from the water Cal The hy- 
draulic diffusivity 7; is defined as Gao} [1987 [Lu] 
eal 2019): 

ki X5 
“a ’ 4 
ae (4) 


where Ay = kyri/ [i saul the fluid seta 


Peres and anes (Barreto Jr. et al.| 
Bela et. al.| ,|Silva et al.| , ponedenne the 


suitable fluid ue flows in region 7; that is, oil or 
water. 

The initial condition, inner and outer boundary 
conditions for the problem are given by: 













Initial Condition (IC): 


Api(r,t = 0) = 0. (5) 
Inner Boundary Condition (IBC): 





lim 2rhkyAy (| = 6(), (6) 
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Figure 6: Reservoir lateral view for rr(t) > rs. 


where 6(t) denotes the Dirac delta function. The 
IBC defined in Eq. (6) represents an instant injec- 
tion of a unitary volume of water. 


Outer Boundary Condition (OBC): 


lint. Apa? = F3,2) = 0: (7) 

Te—oo 
Note that the problem defined by Eqs. to 
remains ill-posed so far, as the PDE presented in 
Eq. demands two spatial conditions per region. 
The missing spatial conditions are obtained from 
pressure continuity and mass conservation at the 
interface between aia ae — condi- 
— pues eecAS or an 


P0ot| Viana) BOD): 


Pressure CCR: 


Ap,(r > r;,t) 
Flow-rate CCR: 


= Apiai(r > 7} ,t). (8 


NS 














where 7 € {1,2}. Eqs. and (9) are required 
to assure that pressure and flow-rate profiles along 
the reservoir are continuous. 

Applying Laplace transform to the problem de- 
fined by Eqs. to (9). the following ODE is ob- 
tained: 





10 (52) — udp; 0, (10) 


r Or Or "i 
for 7 € {1,2,3}, where u denotes the Laplace vari- 
able. 
The general solution of Eq. (10) is given by 


(Lefkovits et al] [1961] [Prats et al.) 2020): 


Ap,(r,u) = alo (+=) +b, Ko (+=) . (11) 


Taking the derivative of Eq. (11) with respect 


to the radius (Abramowitz and Stegun| |1964): 


El Gy)-m (V9 
Or "i ut "i 
(12) 
Relations between the coefficients a; and 6; in 
each region may be obtained from the boundary 
and coupling conditions. The outer boundary con- 
dition defined in Eq. (7) implies that a; = 0 
(Ehlig-Economides and Joseph} 
and Stegun| 


From Eqs. (6) and (12), the inner boundary con- 
dition in Laplace domain may be written as: 


al; (r. =| — by ky (1 =) = 
a a (13) 
= 1 mM 1 
tw Vou QrhkyAz 


From Eq. and the CCRs defined in Eq. (8). 
pressure equality at the interface between regions 
implies that: 













ailo (=) + bio (nf) = 
UD h 
—ai+4ilo (m / ~) — bi41Ko (m / ~) = 0, 
M41 Mi+1 


(14) 


for i € {1, 2}. 
Finally, from Eqs. (11) and (9). mass conserva- 
tion at the interface between regions is expressed 














as: 
kids kids 
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tee een Ky (x 2 ) = 0, 
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(15) 
for i € {1, 2}. 


Writing Eqs. to as a linear system: 





ay - are 
by 0 
[M] |a2} = 0 , (16) 
by 0 
bs 0 


where the first row of matrix M is built using Eq. 
(13), while the second and third rows of matrix M 
are built using Eq. and the fourth and fifth 
rows of matrix M are built using Eq. (15). 

Solving the linear system defined in Eq. (16), the 
coefficients a; and b; in each region may be deter- 
mined. Thereby, pressure change at the wellbore 
can be evaluated as: 


Gi(u) = azlo (roy/=} +b, Ko (rf p(T 


Eq. represents the pressure response mea- 
sured at the wellbore due to the instant injec- 
tion of a unitary volume of water into the reser- 
voir at t = 0, as defined in the IBC given by 
Eq. (6). Thus, the term G,, denotes the im- 
pulse function (or Green’s function) in Laplace do- 
main (Carslaw and Jaeger 
Ramey Jr.||1973), as a consequence of an instanta- 
neous water injection of unitary volume. Note that 


the impulse function at the wellbore must be evalu- 
ated using the properties corresponding to Region 
1, which is the region that contains the well. Thus, 
pressure change during an injectivity test may be 


determined as (Everdingen and Hurst} |1949} 
2020): 


ADw(U) = Ging (U)G@u(u). (18) 
Finally, wellbore storage effects may be easily in- 


corporated into Eq. as follows (Everdingen 
1949| [Ehlie-Economides and Joseph 


Apu s(u,C = 0) 


Ap ne = aa ’ 
Pw ft (U ) 1+ u2CAp,,;(u, C = 0) 


(19) 





where C’ stands for the storage coefficient. 

In multilayer reservoirs, pressure change may be 
determined by combining Eq. (17) with the algo- 
rithm developed by (1991). 

Then, the following procedure may be applied 
to develop a computational implementation that 
computes the pressure change using the formula- 
tion detailed in this Section, for injectivity tests in 
multilayer reservoirs with vertical wells: 


e Define an initial guess for layer flow-rates (e.g., 
according to layer flow capacities) 


e For each time step: 


1. For each layer: 


(a) Compute the water front radius, us- 


ing Eq. 


(b) Define each region of the radially 
composite reservoir model according 
to the water front position 

(c) Build matrix M 

(d) Evaluate the impulse function 


Gyj(u), as defined in Eq. 


2. Compute pressure change in Laplace do- 


main Ap,,- as proposed by |Kuchuk and 
(1991) 

3. Update layer flow-rates in Laplace do- 
main @; as proposed by |Spath et. al. 
(1994) 


4. Apply|Stehfest} (1970) algorithm to invert 


pressure change and flow-rates to the real 
time domain 


e Repeat steps 1 to 4 for the next time step 


2.2. Application to Single-Layer Reservoirs with 
Horizontal Wells 

In reservoirs with horizontal wells, flow in each 

coordinate axis affects the pressure change behav- 

ior. The three-dimensional diffusivity equation for 


an isotropic reservoir is given by (Odeh and Babu 
1990} 2000): 


Ap O?Ap Od?Ap 10Ap 
Om Oy? of n Ot 
There already exist analytical models for single- 
phase flow in single-layer reservoirs with horizon- 
tal wells 
(1989). For a uniform influx wellbore, the 


solution in real time domain is obtained by apply- 
ing (1936) product and impulse functions 


(Gringarten and Ramey Jr.| |1973} |Ozkan et. al.| 





0, (20) 
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(21) 


where Ap,, Ap, and Ap, denote the respective im- 
pulse functions in x-, y- and z-directions. 

When it comes to water injection in single- 
layer reservoirs with horizontal wells, an analyti- 


cal model was developed by |Peres and Reynolds 
(2003), based on the flooding patterns proposed by 


(1961). As explained by|Peres and Reynolds 


, the flood front initially propagates in the 
vertical plane perpendicular to the wellbore axis. 






Therefore, this work proposes to approximate an 
injectivity test in a single-layer reservoir with hor- 
izontal well by a radially heterogeneous reservoir 
model. The reservoir is assumed to be isotropic 
(hie = hy = Ky = A)> Wig. [7] shows a schematics 
of a single-layer reservoir with a horizontal injector 
well, for a given moment when the flood front is 
beyond the damaged region and assuming that the 
wellbore is parallel to the y-axis. Fig. |8| exhibits 
the view in zx-plane of this reservoir model. 

As evidenced in Fig. |7| the approximated radi- 
ally heterogeneous model for injectivity tests with 
horizontal wells considers that the flood front prop- 
agates only along the zx-plane. This assumption, 
albeit unrealistic, is expected to be reasonable, con- 
sidering that the wellbore length is typically much 
longer compared to the reservoir thickness. Be- 
sides, water flow beyond the wellbore tips at early- 
times is also unaccounted for in the model proposed 
by [Peres and Reynolds (2003). 

One should notice, however, that waterfront 
propagation in the horizontal plane becomes in- 
creasingly more relevant as the injection goes on 





Figure 7: Horizontal injector well model (lateral view) 







(2008) (Pores 


i | . Thus, the radially 
heterogeneous model exhibited in Figs. [7] and [8] is 
valid in an approximate sense, considering that the 
flood front has not yet reached a vertical boundary. 

On the other hand, pressure transient front does 
not depend on the flood front position. Hence, 


considering the notation defined by 
Reynolds! (2003), the approximate radially hetero- 


geneous reservoir model shown in Figs. [7] and [8] is 
applicable for the ”first radial - first radial”, ” first 
linear - first radial” and ”second radial - first ra- 
dial” flow-regimes. 

Under those assumptions, and defining r? = 
x? + z’, pressure change during injectivity tests in 
single-layer reservoirs with horizontal wells is gov- 
erned by the following PDE: 





2 : 2 : 
O*Ap; 10 (Se) i OAp; =0. (22) 


Oy? : ror Or nH; Ot 
Then, considering the impulse function defined 
in Section this work proposes that pressure 
change at the wellbore in single-layer reservoirs 
with horizontal wells may be evaluated as: 


t 


Avy, = trust) =4 / Ap,(y, 1) £7 {Gu} (nar, 


0 


(23) 





Figure 8: Horizontal injector well model (view of the za- 
plane) 


where G,,, is the impulse function defined in Eq. 
and Ap,(y,7) stands for the impulse func- 


tion in y-direction, defined as (Daviau et al.||1988| 
Ozkan et al.||1989): 


sivn-for(Et) -a()] 
(24) 


where L indicates the wellbore length. The images 
method can be applied to depict the reservoir verti- 


cal boundaries (Clonts and Ramey Jr.| 
and Carvalho] [1989). 

Eq. (23) consists of an unorthodox application 
of product, since it is derived from 
a non-homogeneous reservoir model, as shown in 
Figs. [7|and[8] A discussion regarding the accuracy 
of Eq. is made in Section A solid proof 
for why Eq. would be an interesting topic for 
further research but is out of the scope of this work. 









3. Results and Discussion 


A comparison was made between a commercial 


flow simulator (CMG| |2010) and the formulations 
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Figure 9: Relative permeability data 


proposed in Section [2|in order to evaluate its accu- 
racy on a set of synthetic cases. Additionally, the 
existing solutions for injectivity tests in multilayer 


reservoirs with vertical wells (Barreto Jr. et al.| 
2011\|Mastbaum et al.|/2021) and single-layer reser- 
voirs with horizontal wells (Peres and Reynolds} 


2003) were also compared to the models presented 
in Section 

In all cases, a piston-like fluid displacement was 
considered. The injection period was set as 96 h. 
Relative permeability data may be seen in Fig. |9} 


3.1. Results for Multilayer Reservoirs with Vertical 
Wells 


The formulation proposed in Section was im- 
plemented and compared to a commercial flow sim- 
ulator (2010). A cylindrical reservoir grid 
was employed, with 70 blocks in the radial direc- 
tion and 1 block in the angular direction. In the 
vertical direction, each grid block represented a dis- 
tinct layer. The reservoir outer radius was set as 
4 km. The innermost grid block is 0.1 diameter 
wide and blocks get coarser as the distance to the 
wellbore increases. 


The solutions developed by |Barreto Jr. et al. 
(2011) and |Mastbaum et al.) (2021) for injectivity 


tests in multilayer reservoirs were also implemented 
and compared to the suggested model. Table 
displays the reservoir properties for the multilayer 
cases. In all cases, a constant injection flow-rate 
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Table 1: Reservoir properties for the multilayer cases 


of 500 m3/d was considered and wellbore storage 
effects were neglected. 


Fig. displays the loglog graphs of pressure 
and pressure derivative for Case A. Pressure data 
are represented in black and derivative curves are 
shown in blue. The formulation suggested in this 
work is indicated by circles. Triangles denote the 
numerical simulator. The dashed lines stand for 


the model proposed by (2021) and 


the solid lines represent the solution developed by 
(2011). The horizontal dashed 
line indicate the theoretical derivative level associ- 
ated with oil properties, while the horizontal line 
in dash-dot pattern stands for the water properties 
theoretical derivative level. 


A great agreement between all methods is ob- 
served. The sharp derivative drop that occurs at 
early times reflects the existence of skin effects 
(Peres and Reynolds 2003 Bela eal, 20208) Siva 
et al.| 2021). Pressure derivative considering the 
numerical simulator data shows some divergence 
from the analytical solutions during this period. 
This may be related to the numerical grid dis- 
cretization. Besides, the piston-like fluid displace- 
ment is more accurately depicted by the analytical 
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Figure 10: Pressure and derivative curves for case A 


models than the flow simulator. Still, after the 
damaged region is swept by water, all derivative 
curves converge. 

It is also important to highlight that all meth- 
ods attained the theoretical derivative level associ- 
ated to the water properties. This indicates that 
the proposed formulation may be applied to esti- 
mate the reservoir equivalent permeability, based 
on the constant derivative level (Banerjee et al.| 
[2020b). Fig. |10| also shows 
that the derivative curves attain values close to the 
oil derivative level at early-times but quickly de- 
taches from this theoretical level, as the derivative 
drop starts. 

Fig. presents the layer flow-rate profiles for 
each methods, in semilog scale. The black curves 
denote layer 1, while blue curves stand for layer 2 
and red curves represent layer 3. Fig. [I1]evidences 
that layer flow-rates at the wellbore change in time. 
Hence, as mentioned in Section[I] it is unfeasible to 
assure that flow-rate within the swept zone remains 
constant in multilayer reservoirs. Layer 1 presents 
the lowest flow capacity, while layer 3 presents the 
highest flow capacity. Thereby, flow-rate in layer 
1 is lower than in layer 2, which is lower than in 
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Figure 11: Layer flow-rate profiles for Case A (layer 1 in 
black; layer 2 in blue; layer 3 in red) 


layer 3. 

As disclosed in Fig. the proposed formu- 
lation exhibits a close convergence to the simu- 
lated data and the analytical model from 


(2021). However, flow-rate profiles 


computed using the analytical solution developed 
by|Barreto Jr. et al.] diverged from the other 
methods in layers 1 and 3. This result is partic- 
ularly interesting, since accurate knowledge over 
layer flow-rate data is crucial for parameter estima- 


tion techniques (Coutinho et al.| /2010) 
2077 Bela ot al BU20E). 

Nonetheless, it is interesting to notice that the 
formulation proposed by 
was able to provide pressure and pressure deriva- 
tive curves that accurately match the other meth- 
ods, despite the divergence observed in the layer 
flow-rate profiles. 

Pressure and derivative curves for case B may be 
seen in Fig. Case B presents the same reser- 
voir properties as Case A, expect for oil viscosity. 
Therefore, the theoretical oil derivative level is now 
lower than the water derivative level. The hump 
in pressure derivative curves that occurs between 
t = 0.01 h and t © 0.2 h reflects the existence of 
skin effects in reservoirs where flow is favorable to 
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Figure 12: Pressure and derivative curves for case B 


11 


oil (Bela ct al,] 2020n] Silva ct al] 2021). After this 


derivative hump, flood front propagates through- 
out the undamaged zone. 

Similarly to Case A, all methods presented a 
close agreement. The three analytical formula- 
tions converged during the entire test, while numer- 
ical data presented some slight divergences at early 
times, when the flood front is still within the dam- 
aged zone. Despite that, the overall agreement is 
quite good and all methods converged towards the 
theoretical water derivative level, indicating that 
late-time data can be used to determine the reser- 
voir equivalent permeability. 

Layer flow-rates for Case B are reported in Fig. 
The formulation suggested in Section[2]exhibits 
an excellent agreement with the solution developed 
by (2021). Flow-rates evaluated 
using these analytical models show a similar be- 
havior to the numerical simulator flow-rate curves. 
Yet, in layers 1 and 3, the match between analyti- 
cal and simulated data is not as close as observed 
in Case A. 

On the other hand, flow-rate profiles determined 


by the solution proposed by |Barreto Jr. et al. 
(2011) presented once again noticeable differences 
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Figure 13: Layer flow-rate profiles for Case B (layer 1 in 
black; layer 2 in blue; layer 3 in red) 


compared to the other methods. Figs. and 
show that the solution from |Barreto Jr. et_al. 
provided a good match for pressure data but 
yielded layer flow-rate profiles that diverged from 
the other models, as occurred in Case A. 


Pressure and derivative data for Case C may be 
observed in Fig. The overall behavior is sim- 
ilar to Case A: pressure change presented a close 
agreement for all methods, as well as the deriva- 
tive curves. The blunt derivative drop at early- 
times indicates the existence of formation damage. 
At late-times, the theoretical derivative level as- 
sociated with water properties is attained by all 
models. 


In Case C, layer flow-rate curves, portrayed in 
Fig. |15} provide a more interesting analysis than 
pressure and derivative data. Fig. [15]evidence that 
both the formulation proposed in Section and 
the solution presented by 
exhibited a great agreement with the numerical 
data. As mentioned earlier, accurate evaluation 
of flow-rates in each layer is helpful to some in- 
terpretation techniques that aim to estimate in- 


dividual layer properties (Ehlig-Economides and 
1987| [Coutinho ct al. [2010] 
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Figure 14: Pressure and derivative curves for case C 
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(Guimaries] BOTT) [Beta et al} 20200). 


Similar to the previous cases, flow-rates com- 
puted using the solution from [Barreto Jr. et al.| 
diverged from the other analytical models 
and from the numerical simulator. Moreover, the 
solution from presented 
the most significant variations in layer flow-rates, 
even though this is the only analytical model that 
requires flow-rate to be constant within the swept 
region. Despite that, it is interesting to notice 
the close agreement presented by the pressure and 
derivative curves in Fig. 

Pressure and layer flow-rate profiles displayed in 
Figs. to [15] were employed to estimate layer 
permeabilities using the Delta Transient Method 
(Galvio_and _Guimaraes| 
Galvao} |2017). This procedure requires pressure 
and flow-rate data from only two distinct time 
steps, which is a significant operational advantage 
0200) 

In this work, the late-time logarithmic approx- 
imation based on water properties was employed 
to evaluate layer permeabilities. For all methods, 
pressure and flow-rate data at t = 6.06 h and 
t = 60.6 h were used to evaluate layer permeabil- 
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Figure 15: Layer flow-rate profiles for case C (layer 1 in 
black; layer 2 in blue) 


























This Work Num. Simulator | Mastbaum et al. | Barreto et al. 

Case True | Est. Er. (%) | Est. Er. (%) | Est. Er. (%) | Est. Er. (%) 
500 | 498 -0.3 493 -1.4 498 -0.3 429 -14.2 
A 600 | 599 -0.2 589 -1.8 5099 -0.2 560 -6.7 
700 | 700 0.0 691 -1.3 700 0.0 767 9.6 
900 | 501 0.1 500 0.1 001 0.1 481 -3.7 
B 600 | 600 0.1 5099 -0.2 600 0.1 590 -1.6 
700 | 700 -0.1 700 0.0 700 -0.1 716 2.2 

C 500 | 499 0.2 497 -0.7 499 -0.2 448 -10.4 
700 | 701 0.1 683 -2.4 701 0.1 804 14.8 














Table 2: Estimated Layer Permeabilities (Values in mD) 


ities. Figs. and evidence that, for all 


cases, a radial flow regime occurs during this time 
range, since derivative remains constant. 

Table [2] displays the estimated layer permeabil- 
ities considering each approach employed to com- 
pute the pressure change. A great accuracy may 
be observed in all cases. Moreover, the formulation 
presented in this work yielded values very close to 
the estimates obtained from the flow simulator and 
the model proposed by (2021). 
The solution developed by (2011), 


in its turn, showed some differences compared to 
the other models (particularly in case B), despite 
presenting a decent accuracy. 

These divergences are likely related to the mis- 
match between layer flow-rate data. As displayed 


in Figs. and flow-rates evaluated using 


the solution developed by |Barreto Jr. et al.) (2011) 


diverged from the other approaches, which yielded 
similar flow-rate profiles. These results combined 
with the higher errors in layer permeabilities ob- 
served in Table |2|suggest that layer flow-rate pro- 
files computed using their model are less accurate 
than the flow-rate curves obtained from the other 
methods. 


3.2. Results for Single-Layer Reservoirs with Hor- 
izontal Wells 

A computational implementation of the analyt- 

ical model presented in Section [2.2] was developed 
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and compared to the solution proposed by 
for injectivity tests in single- 
layer reservoirs with horizontal wells. Results were 
also compared to the commercial flow simulator 
B70). 

A cartesian grid composed of 61 x 61 x 9 blocks 
was employed. Again, blocks closer to the well- 
bore were more refined than blocks further from it. 
Reservoir grid was 6,000 m wide (x-direction) and 
6,300 m long (y-direction). The wellbore was built 
parallel to the y-direction. Length (y-direction) 
of each block containing the wellbore was set as 
100 m, while width (z-direction) and height (z- 
direction) were set as 0.89 m. Therefore, the skin 
zone radius was set as 0.5 m in the analytical mod- 
els, so that the damaged region presents the same 
cross-sectional area that the numerical simulation 
grid. A hybrid refinement of 4 (radial direction) by 
4 (angular direction) by 4 (well axis direction) was 
employed at all grid blocks containing the wellbore. 

The numerical simulator considers an infinite 
conductivity wellbore, while both analytical solu- 
tions assume uniform influx. Either the average 


pressure technique (Kuchuk et al. |1991} 
(Chen et. al.|}1991) or the equiv- 
alent pressure point (Ozkan et al.| |1989| 
may be applied to obtain the 


pressure change at an infinite conductivity well- 
bore using a formulation that assumes uniform in- 
















ithe k h L kg Ts 
Case 

(cP) (mD) (m) (m) % (mD) (m) 
D 48 2000 25 700 0.25 £500 0.5 
E 1.1 2000 25 700 O.25 £500 0.5 








Table 3: Reservoir properties for cases with horizontal well 


flux. While the average pressure technique may 
be more accurate, its implementation is also more 
complex. For this reason, the equivalent pressure 
point was employed to obtain the pressure change 
corresponding to an infinite conductivity wellbore. 
Thus, pressure change in the analytical formula- 
tions was evaluated at x = z = ry and y = 0.732 
times the wellbore length, as suggested by 
(1080), 

Layer properties for the tested cases may be 
found in Table 3} Injection flow-rate was defined as 
5,000 m3/d for Cases D and E. As in Section [3.1] 
wellbore storage effects were disregarded. Reser- 
voir properties are the same in both cases, except 
for the oil viscosity. The horizontal well was as- 
sumed to be at the center of the layer; that is, no 
off-centering was considered. Waterfront radius at 
the end of the 96h injection period, determined by 
Eq. (1). was equal to 8.8 m. Thus, the flood front 
has not yet reached the reservoir vertical bound- 
aries and the formulation developed in Section 
is applicable. 

Fig. displays the pressure and pressure 
derivative data for Case D. Black curves denote 
the pressure data while the blue curves correspond 
to the derivative profile. The model proposed in 
Section is indicated by circles. Triangles rep- 
resent the numerical simulator and the solid lines 
correspond to the solution presented by {Peres and| 
(2003). The horizontal dashed line in- 
dicates the theoretical pressure derivative profile 
for ’first radial - first radial” flow regime, as de- 
fined by [Peres and Reynolds] (2003). At late-times, 
the flood front is still propagating in the vertical 
plane, while pressure diffusion is occurring in the 
horizontal direction. Therefore, a ”second radial 
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Figure 16: Pressure and derivative curves for case D 


- first radial” occurs at the end of the test. The 
dash-dotted horizontal line stands for the theoret- 
ical level for this flow regime. 

A good agreement is observed between all curves. 
As in Section [3.1] the pressure derivative drop ob- 
served at early-times indicates the existence of for- 
mation damage. Pressure derivative for the numer- 
ical simulator presents some oscillations for short 
times (¢ < 3-107? h). These oscillations are 
possibly explained by the grid discretization and 
how the numerical simulator accuracy to represent 
a piston-like flow. Despite that, numerical data 
match well the analytical models. 

It is worthy to mention that, interestingly 
enough, the proposed formulation based on an 
unusual application of product 
presents a closer match to the simulator derivative 
curve than the solution developed by {Peres and] 
at early-times. Furthermore, all 
methods attain the late-time theoretical derivative 
level. 

The results for Case E may be seen in Fig. 
All methods exhibit a great match not only for 
pressure but also for derivative data. In Case E, 
the derivative curves quickly detach from the early- 
time theoretical level. Still, the derivative signa- 


Pressure and Pressure Derivative Profile 
1 1 1 : 





10° 





© This Work 
V_ Flow Simulator 
Peres and Reynolds 





























107 
10° 


10° 10! 10° 


Figure 17: Pressure and derivative curves for case E 


ture that indicates a late-time radial flow is clearly 
observed. This constant derivative level is useful 
to determine the reservoir permeability (Odeh and 


1990} |Rosa and Carvalho} |1989| 
and Reynolds} |2009)). 


Since Case E presents a mobility ratio lower than 
unity, no derivative drop is observed. One could 
expect a derivative hump instead, similar to cases 
B and C. However, this hump may be masked by 
the derivative rise due to the linear flow regime. 
Yet, it is interesting to notice once again the close 
agreement presented by the proposed solution, con- 


sidering the application of (1936) product 


for a non-homogeneous reservoir model. 


4. Summary and Conclusions 


This work presents an alternative solution for 
injectivity tests in single-layer reservoirs. The 
suggested formulation uses a radially composite 
reservoir model with moving boundary to depict 
the swept region (Bratvold and Horne [1990] 
and applies impulse functions to obtain the pres- 


sure change at the wellbore. The single-layer solu- 
tion may be coupled to the algorithm developed by 
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Kuchuk and Wilkinson (1991) to evaluate the pres- 


sure behavior in multilayer reservoirs. A formula- 
tion is also proposed for injectivity tests in single- 
layer reservoirs with horizontal wells by combining 
the proposed single-layer solution to an unusual ap- 


plication of (1936) product. 


The suggested formulations were compared to 


a commercial flow simulator (CMG} |2010) and to 
previously existing solutions 
(2021). Results for the multilayer cases (seen in 
Section showed that the proposed model fits 
well the pressure and flow-rate data provided by 
the numerical simulator. 

Furthermore, results evidence that layer-flow 
rates do not remain constant in time. This is 
a relevant outcome, as it demonstrates that, rig- 
orously speaking, the 
steady-state theory does not hold for mul- 
tilayer reservoirs. Despite that, it is interesting to 
notice that all methods presented a good match 
for pressure and pressure derivative data. On the 
other hand, layer flow-rate curves computed by the 
solution developed by|Barreto Jr. et al.] pre- 
sented significant divergences compared to the an- 
alytical model proposed in this work, to the other 
models. 


The Delta Transient Method (Galvao and 
2017) was employed to estimate indi- 


vidual layer permeabilities. Results showed a satis- 
factory accuracy. Furthermore, estimates obtained 
from all models except the model presented by [Bar-| 
were very close to each other. 
This outcome and the aforementioned divergences 
in layer flow-rate curves suggest that the formula- 
tion developed by is less 
accurate to evaluate flow-rate behavior than the 
other approaches, despite yielding a decent match 
for pressure and derivative data. 

Results from Section |3.2|verify that the unortho- 
dox application of product sug- 
gested in Section [2.2| was able to provide a formu- 
lation for injectivity tests in single-layer reservoirs 


with horizontal wells that matches well the results 
from a commercial flow simulator and also the so- 


lution developed by |Peres and Reynolds} (2003). 


A solid proof to justify the application of 
product using a radially heterogeneous 
reservoir model was out of the scope of this work, 
and remains as suggestion for further research. 


Nomenclature 


a = Multiplicative coefficient, defined in Eq. (11) 
b = Multiplicative coefficient, defined in Eq. (11) 
C = Wellbore storage coefficient [m?/Pa, STB/psi, 
m®/ (kgf /em?)] 

c = Compressibility, [(Pa)~!, (kgf/cm?)~"] 

erf = Error function 

fi, = Fractional flow derivative [-| 

Gy = Impulse function in Laplace domain 

h = Thickness |m, ft] 

I; = Modified Bessel function of first kind and 
order 7 

K; = Modified Bessel function of second kind and 
order 7 

k = Permeability [m?, mD] 

k, = Relative permeability [-] 

L = Wellbore length |m, ft] 

M = Matrix defined in Eq. 

p = Pressure [Pa, psia, kef/cm?] 

q = Flow-rate [m?/s, STB/d, m?/d] 

r = Radius [m, ft] 

Tw = Wellbore radius [m, in] 

Sw = Water saturation |-] 

t = Time |s, hj 

u = Laplace variable 

6 = Dirac delta function 

n = Hydraulic diffusivity coefficient [m?/s, mD- 
psia/cP, mD- kgf/cm?/cP] 

= Fluid mobility [(Pa-s)~1, (cP)~"] 

@ = Porosity [-] 

T = Silent time integration variable |s, h] 
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Subscripts 


F = Related to the flood front 

i = Related to region i 

inj = Related to the injected fluid 
j = Related to layer j 

o = Related to oil 

s = Related to the damaged region 
i= Total 

x = With respect to x-direction 

y = With respect to y-direction 

z = With respect to z-direction 
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